Coupled channel analysis of the p meson decay in lattice QCD 
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We employ a variational basis with a number of qq and tttv lattice interpolating fields with quan- 
tum numbers of the p resonance to extract the discrete energy spectrum in a finite volume. In the 
elastic region, this spectrum is related to the phase-shift of the continuum scattering amplitude by 
Liischer's formula and the relation allows the extraction of resonance parameters from the spectrum 
calculation. The simulations are performed at three different total momenta of the coupled qq — nn 
system, which allows us to extract the p-wave scattering phase at five values of pion relative mo- 
menta near the resonance region. The effective range formula describes the phase-shift dependence 
nicely and we extract the resonance mass m p — 792(7) (8) MeV and the coupling g p7I -ir = 5.13(20) at 
our m„ ~ 266 MeV. The coupling ,g pw ^ is directly related to the width of the p meson and our value 
is close to the value derived from the experimental width. The simulations are performed using dy- 
namical gauge configurations with two mass-degenerate flavors of tree-level improved clover-Wilson 
fermions. Correlation functions are calculated using the recently proposed distillation method with 
Laplacian Heaviside (LapH) smearing of quarks, which enables flexible calculations, in many cases 
with unprecedented accuracy. 

PACS numbers: 11.15.Ha, 12.38.Gc 



I. MOTIVATION AND INTRODUCTION 

Almost all hadrons listed in the Particle Data Group 
PJ tables are unstable, most of them decaying strongly. 
In quenched calculations, where vacuum quark loops are 
disregarded, all hadronic states appear as stable states. 
In full QCD, on the other hand, truly asymptotic expo- 
nential behavior is always dominated by the lowest stable 
end product. This is unsatisfactory. 

In continuum physics experiments resonances are iden- 
tified via the scattering cross section and subsequent 
phase-shift analyses. In the lattice discretization of QCD, 
instead, one studies the correlation functions of hadron 
interpolators for Euclidean time distances. The result is a 
combination of exponentially decaying terms, each corre- 
sponding to the energy level of a contributing eigenstate. 
Due to the finiteness of the lattice system, the energy 
levels are discrete. The spectral density is related to a 
discretization of the cross section. However, in realistic 
lattice simulations only very few such levels can be de- 
termined. The typical gaps are 0(2ir/L) for lattices of 
spatial extent L; for most simulations this corresponds 
to level spacing 0(400) MeV. 

However, as has been pointed out in a seminal paper 
by Luscher 0, HJ, for a resonating system the discrete 
spectrum obtained in a finite volume can be related to 
the phase-shift of the continuum scattering amplitude in 
the elastic region. The resulting volume dependence of 
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the spectrum can then be used to explore the resonance 
properties Q ■ Model simulations in two dimensions Q as 
well as in four dimensions @ demonstrated the feasibility 
of that approach. The original derivation in the decaying 
particles rest frame was then extended to moving frames 
[7H9( , thus enhancing the practical applicability, allowing 
one to obtain the phase-shift at more momentum points 
for a given lattice size. Because of several problems there 
have only been a few attempts to apply that scheme to 
the decay p — > rnr P4l3| , while the first lattice estimate 
of the p — > 7T7T amplitude [l5[ did not apply Liischer's 
method. Note that widths for most of the other reso- 
nances have not been determined on the lattice at all. 

There are two major complications. The first one con- 
cerns the hadronic lattice interpolators used. Let us 
assume that we work with the fully dynamic vacuum, 
i.e., including the dynamical quark vacuum loops in a 
full QCD simulation. Naively one would expect that, 
even if one correlates only quark- antiquark interpolators 
with the correct quantum numbers of the p, due to the 
vacuum loops, tttt intermediate states should also con- 
tribute and affect the energy levels accordingly. This 
is hardly observed; actually already in model calcula- 
tions Q it proved necessary to include both, the heavy 
boson and the two light bosons in the set of interpola- 
tors. Similar observations were made in other calcula- 
tions involving baryon and meson correlation functions 
[l6l - [l9| . The obvious interpretation is that the overlap of 
the quark-antiquark interpolators with the meson-meson 
decay channel interpolators is too weak to have been ob- 
served. 

For that reason one should extend the set of hadron 
interpolators to include both, various versions of the 



2 



quark-antiquark interpolator (like, e.g., different Dirac 
structure or different quark smearing functions), as well 
as meson-meson interpolators. The latter involve four 
propagating fermions and the corresponding entries of 
the correlation function usually will involve backtracking 
loops. In addition to this technical complication there 
is also the notorious issue of statistical weight for such 
contributions. The so-called distillation (or Laplacian- 
Heaviside quark smearing) method introduced in [20[ and 
employed in [l9|, [2ll - [23 | helps us significantly to deal with 
that problem. 

The second challenge concerns the energy levels. One 
works with several hadronic interpolators, all with the 
correct quantum numbers and total momentum in the 
given channel. The diagonalization of the correlation 
matrix gives the eigenstates and eigenenergies according 
to the so-called variational method |25l - [28| . The set of 
lattice interpolators should be large enough to be able to 
represent the leading eigenstates and thus the leading en- 
ergy levels. The better the set is, the better the results 
will be and the more energy levels can be determined, 
depending of course also on the available statistics. In 
previous calculations aimed at p meson decay, at most 
two interpolators were used: one quark-antiquark and 
one pion-pion interpolator. We extend this to a larger 
interpolator basis. 

For our calculation we use one lattice ensemble with 
rif = 2 dynamical mass-degenerate light quarks and 
clover-improved Wilson fermionic action (generated in 
context of the work [2^, [3(| in order to study reweighting 
techniques). The ensemble consists of 16 3 x 32 lattices 
with spatial extent 1.98 fm and m T ~ 266 MeV. We con- 
sider cross-correlations of several interpolators (16 for the 
p channel, 6 for the pion channel) and solve the gener- 
alized eigenvalue problem to reliably determine the two 
lowest energy levels. We study the p channel for three val- 
ues of the total momentum and obtain the elastic phase- 
shift in the resonance region. 

Section HQ gives an overview of the methods: quarks 
sources, interpolators, variational analysis, phase-shift 
relations and finite time effects. In Sect. IIIII the set of 
configurations and details on the computations are sum- 
marized and in Sect. IIVI we discuss the results: correla- 
tion functions, energy levels, phase-shift and resonance 
parameters. 

Reference [3lJ suggests an alternative approach which 
has recently been investigated in (32j. Furthermore an- 
other procedure has been suggested in [33| . 



II. TOOLS 

A. Phase-shift formulas, brief review 

On finite lattices there are, strictly speaking, no 
asymptotically free states and the energy spectrum is al- 
ways discrete. It was pointed out by Luscher @, 0] that, 
assuming a localized interaction range, the energy level 



of a correlation matrix for channels with resonances in a 
finite volume can be related to the corresponding phase- 
shift in infinite volume in the elastic region (i.e., where 
only one decay channel is open). The relation was de- 
rived for interpolators with spatial momentum zero. For 
a particle like the p meson, which can decay into two pi- 
ons with back-to-back momenta, the available momenta 
are discrete on finite lattices and depend on the spatial 
extent. 

In the noninteracting case the various two-pion energy 
levels will decrease with growing volume and this leads 
to level crossing with the stable p state. If interaction is 
switched on the level crossing is avoided and the energy 
levels "change their identity" . This was demonstrated in 
a two dimensional resonance model in J5| as well as in 
four dimensional </> 4 -model simulations [6j. 

For the analysis of resonances in that method one needs 
several ingredients. The set of interpolators should over- 
lap with both, the single particle content (i.e., for a me- 
son the quark-antiquark component) as well as the two 
particle content (i.e., the meson- meson decay channel). 
Furthermore it should be possible to analyze more levels 
than just the ground state energy. Third, in the origi- 
nally proposed method one needs several spatial volumes 
to obtain the phase shift at several values of relative mo- 
mentum. This makes the approach costly. 

The third aspect can be ameliorated, though, by study- 
ing also channels with nonvanishing total momentum 

P = — d with deZ 3 . (1) 

Lj 

In our simulation we study the cases 

d= (0,0,0), (0,0,1), (1,1,0) (2) 

and permutations, which have previously been combined 
in the simulation [l3[. Different values of P allow to 
obtain the phase shifts at different values of pion relative 
momenta. The lowest ttt: state in the p channel with 
|P| = is 7r(27r/L) 7r(-2vr/L) (due to I = 1) and is 
significantly above the p resonance in typical simulations. 
In the case of a p with |P| = 2-k/ L, the 7r(0) ir(2ir/L) is 
closer to the resonance region, for example. However, 
this case involves relativistic kinematics in the nonzero 
momentum frame as pointed out in @-|9| • The relativistic 
distortion reduces the full cubic symmetry Oh to that 
of prismatic dihedral groups, i.e., to the symmetry of a 
cuboid (quadratic prism) D^h for total momenta of type 
(0,0,1) and to the symmetry of a rhombic prism £>2/i for 
momentum (1,1,0). 

In the laboratory frame, the total 3-momentum of two 
noninteracting bosons in a cubic lattice of volume I? and 
periodic boundary conditions is 

P = Pi + P2 = ^ ( 3 ) 
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and the energy is 



E = E 1 +E 2 = \Jm 2 + p( + \jm 2 + p 2 
with pi = — rii , rii G Z . 



(4) 



The velocity v = ~P/E gives the relativistic boost fac- 
tor 7 = — v 2 . In the center-of-momentum frame 
(CMF) the total momentum vanishes and the bosons mo- 
menta are 



Pi = -P2 = P ■ 

The energy in the CMF is 

E CM = 2yjm 2 + p* 2 = Eh 



(5) 



(6) 



and the momentum is related to the laboratory frame 
through 



* 1 -if ^ 

P = ^7op (Pi - P2j , 



where the boost factor acts in direction of v, 



(7) 



l op 1 P = Pll/7 + P-L , P|| = v(p- v)/|v| 2 , P±=p-P||. 

(8) 

The relativistic 4-momentum squared is invariant, thus 
the relation to the laboratory energy E is 



E 2 CM = E 2 -P 2 



™* 2 — 1 ™ 2 
P - 1 h CM ~ 



(9) 



the scattering phase-shift. Expressed through the CMF 
variable 



»*2 



2n 



(12) 



one obtains relations of the form tan S(q) = f(q) for tran- 
scendental functions f(q). 

We concentrate on the decay p — > nn where the two 
pions are in p-wave (£ = 1). Details have been discussed 
in the original papers 0, 0, 0-0 EH . For completeness we 
summarize here only the relevant final expressions, where 
phase-shifts are expressed in terms of the generalized zeta 
function defined by 



7 A (*-n 2 \ - V -^*»( x ) 
^■£ml s ><7 I — ( x 2 _ q 2y 



(13) 



x£P d 



P d = jx e K 3 I x = 7" 1 + I) , m G Z 3 J , 

y em ( x ) = |x| £ y £m (x) , 

and y^m are the harmonic polynomials to the spherical 
harmonics functions Y(, m . The zeta function has to be 
analytically continued to s = 1. The simpler form for 
d = is given in Q . A rapidly convergent expression for 
nonvanishing d is derived in (9|. We numerically com- 
pared the different representations of the zeta functions 
of Q and 0] and found agreement. 

The symmetry groups of the sum appearing in 
Zim (jTB")) are Oh, D^h and D2/1 respectively for 
d = (0,0,0), (0,0,1) and (1,1,0). The J p = 1" states 
appear in the specific representations of these symmetry 
groups and the final expressions for the phase-shifts are: 



Due to the coarseness of the lattice we replace in our cal- 
culations this continuum dispersion relation by the lattice 
dispersions relation as suggested in 0, i.e., 



cosh E, 



CAia 



k=l 



cosh^a-2^sin 2 ( ^ ) , (10) 



2sin^-j =2 cosh — C ^ /a — 2 cosh too . (11) 



Zero momentum P = (0, 0, 0) 

(for irrcp Tf in 0,0 0: 



tan 5(q) 



^' 2 q 

2 00 (l;g 2 ) 



Nonzero momentum P = (0,0,1)^ 
(for irrep A 2 in D^h) [7j : 



(14) 



For the interacting case, the momenta of individ- 
ual pions in the laboratory frame are no longer multiples 
of 2tt/L. Assuming a localized interaction region one as- 
sociates the outside region with that of two free bosons. 
The observed energy levels E n are shifted and related to 



tan 6(q) 



77r 3,/2 5 3 



7% d o(i;9 2 ) 



2 2 d (l;9 2 ) 



(15) 



Nonzero momentum P = (1, 1, O)- 2 ^ 

(for irrep B± in D 2 h) [H|: 



tan S(q) = = 7^ 3/ V _ (lg) 

q 2 Z« (l;q 2 ) J\ Zf (l;q 2 ) + uU: (Zf 2 (l;q 2 ) Z«Jl-q 2 )) 
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We independently derived this relation and we agree with 



B. Variational analysis 

To extract the lowest two energy levels with the quan- 
tum numbers I G (J PC ) = 1 + (1 ) of the p meson as 
well as the ground state energies with quantum num- 
bers I G (J PC ) = 1~(0 '") of the pion, we construct 
a matrix C(t)ij of lattice interpolating fields contain- 
ing both quark-antiquark and meson-meson (in our case 
pion-pion) interpolators 

c(t) l0 = j2 e ~ tEn ( \ o >\ n )( n \°»- ( 17 ) 

n 

For this matrix, the generalized eigenvalue problem 

C{t)i^ = A (n) (*)C(io)# n) (18) 

is solved for each time slice. For the eigenvalues \^ n > (t) 
one obtains 

AW(f) oc C ~ tE " (l + O (e- 4 A£ ")) , (19) 

so that each eigenvalue is dominated by a single energy 
at large time separations. This method is called the vari- 
ational method (25T4281 ] . For a detailed discussion of the 
energy difference A.E n , which is in general given by the 
difference between the energy level in consideration and 
the closest neighboring level, please refer to [28j . 

We calculate the eigenvector components of the regular 
eigenvector problem 

C{t Q )-?C{t)C(t )-^ in) ' = A< n >(t)$W . (20) 

In addition to the eigenvalues, the eigenvectors provide 
useful information and can serve as a fingerprint for a 
given state. To track the eigenvalue corresponding to a 
given energy over the full range of time separations, the 
eigenvalues have to be sorted, either by their magnitude 



this expression, originally presented in [j| [l3| . 
I 

or by scalar products of their eigenvectors. In the pres- 
ence of backwards running contributions caused by the 
finite time extent of the lattice, a combination of both 
methods works well: the eigenvalues are sorted by mag- 
nitude at low time separations and by scalar products at 
larger time separation. For our analysis we choose this 
method. 



C. Interpolators 

For the p channel we employ fifteen quark-antiquark 
interpolators and one pion-pion interpolator with J PC — 
1 and \I, I3) = |1, 0) in the variational basis for each of 
the three choices for P as given in ([2]). All previous sim- 
ulations aimed at determining the p meson width used 
at most one quark-antiquark and one pion-pion inter- 
polator and extracted the two lowest energy levels from 
a 2 x 2 variational basis. This may not be reliable if 
the third energy level is nearby and does not allow test- 
ing whether the resulting two levels are robust against 
the choice of interpolators. A larger basis enables us 
to exploit the dependence of the extracted energies on 
the choice of the interpolators. It also indicates whether 
the lowest two states can be reliably extracted using our 
quark-antiquark interpolators alone, or whether the pion- 
pion interpolators are required in the variational basis. 

The 15 different quark-antiquark interpolators Of ype 
(type = 1,..,5, s — n,m,w) differ in type (Dirac and 
color structure) and width of the smeared quarks q s . We 
use three different smearing widths s = n, m, w (narrow, 
middle, wide) for individual quarks and all quarks in a 
given interpolator have the same width s in this simula- 
tion. (Choosing different quark widths within an inter- 
polator is a straightforward generalization and one just 
needs to pay attention that the resulting C-parity is cor- 
rect.) The details on the smearing are given in Subsect. 
Ill El The interpolator 0% is the 7T7t interpolator whose 
structure is explained at the end of this subsection. Our 
sixteen p interpolators are: 
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Ol(t) = 7f u s (x) A tll c iPx u s {x) ~ {u s o d s } (s 



n, m, w) , 



°2W = 72 Us ^ Tt^i7i e iPx u s (x) - <-» 4} (s = 



n,m,w) 



°3 (*) = X! 72 "'W^J ° iPX - {li s O 4} (s = n, m, 10) 



OK*) - E 72 ^ ^< - V 4 C iPx ] Us (.T) - {«. ^ d s } ( 



s = n,m,w) 



0%{t) = X] 71 e ^ A *7j75 5[c iPx ^/ - V/c iPx ]u s (x) - {u s o 4} (s = n,m,w) 



x,i,j,k 



Or n (t) = ^+( Pi )tt-( P2 ) - 7r-( Pl )7r+(p 2 )] , t^Cpi) = ^ ^hs^e*^*) . 



(21) 



In the pion interpolator t^- denote the corresponding 
combination of Pauli matrices and the tttt interpolator 
Og is always composed from narrow quarks. The covari- 
ant derivative (often denoted by 



^(x,y) = £/ i (x ) 0)5 x+ i,y - Uj (x- i, 0)* 



(22) 



is used in some of the quark-antiquark interpolators (used 
already in a number of lattice simulations, e.g. [lfl H3]) 
and will also be employed to prepare smeared quarks q s 
below. It acts on the spatial and color indices and leaves 
time and Dirac indices intact. The linear combinations 
in 04,5 are required for good C-parity. The polarization 
vector A of the quark-antiquark vector current depends 
on the total momentum P 



as 



d= (0,0,0) 
d= (0,0,1) 
d= (1,1,0) 



A = (0,0,1), Pl = -fA, p 2 = ^A. 
A = d, Pi=0, p 2 = P 

A = d, Pi=0, p 2 = P 

(23) 



Our choices for tttt interpolators Oq (|21[) with momentum 
projections for individual pions (|23p are the same as in 

• For d = (0, 0, 0) with the symmetry group Oh 
our interpolator transforms according to the three- 
dimensional representation Tj~ (so just like e z ) un- 
der elements of Oh- 

• For d = (0, 0, 1) with the symmetry group 
Dih the interpolator transforms according to one- 
dimensional A2 (like e z ) under elements of D^h- 

• For d = (1,1,0) with the symmetry group 
D-2h our interpolator transforms according to one- 
dimensional B± (like e x + e y ) under elements of 
D 2 h- Note that the interpolator C 6 with pi = 
(1,0,0) and P2 = (0,1,0) has the same total mo- 
mentum, but it has positive parity and it will not 



appear as an eigenstate for interpolators with B 1 
transformation properties. 

For the isovector pion J PC = h correlation matrix 
we use altogether 6 interpolators, using three smearing 
widths for each of the two Dirac structures, 

X 

Ti=75, r 2 = 757t ; s = n,m,w. (24) 

D. Correlators and contractions 

In the p channel we compute 16 x 16 correlation ma- 
trices for 

C ]k (t f ,U) = (0|O i{t f )Ol{ti)\0) , J, k = 1, 16 , (25) 
where the indices j and k stand for the combination 



(type,s) in Ot VP e,s (EH)- These correlators involve (cf., 
Fig. [T]) connected contractions (a,b) , singly disconnected 
contractions (c), and contractions (d,e). Due to the mo- 
mentum projections at the sink time slices tf, the con- 
tractions (c) and (d) in particular require the propaga- 
tors M _1 from any spatial point at the sink time slice 
t f = 1 V,. 



E. Laplacian Heaviside smearing for quarks and 
the distillation method 



Since calculating all elements of M _1 for the fermion 
Dirac operator matrix M is prohibitively time consum- 
ing, we apply the distillation method proposed in (2ol [ . 
This method is based on a special kind of smearing 
for quarks, that allows treatment of all necessary con- 
tractions. All quarks are smeared according to a pre- 
scription similar to the conventional one q~ auss (x, t) = 
c er 3 v <j,( Xj £) where V 2 denotes the 3D lattice Laplacian 
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(a) 



(b) 



FIG. 1. Contractions for our correlators with qq and 7T7t 
interpolators. 
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acting in a time slice. The major simplification is due to 
the spectral decomposition^ 



N 



./•(A) = ]T/(A( fe >)« 



(26) 



fe=i 



for matrix A = V 2 giving e CT = v2 = £)f e CT * A< *V fe M fe )t . 
Here A^ fc ^ and are eigenvalues and eigenvectors of V 2 
(|22p which is a N^N C x N^N C matrix on a given gauge 
configuration 

VL, xV (*)^W = A (fc) (t)4c ) (*) (27) 

and all the resulting eigenvalues are negative. The choice 
of smearing is arbitrary and instead of this Gaussian 
smearing we use the truncated spectral representation 
of the unit operator (also called the Laplacian Hcaviside 
(LapH) smearing) , as proposed in [2(| and employed also 

in 



N c N'f 



9s ^e(a 2 +V 2 ),= Y, ©(- S 2 + A (fe) ) 



<1 ■ 



k=l 



gT(x,t) = $>&>(*) v$>J(t) q ac '(x',t) 



k=l 

a,a' = l,..,N d =4, c,c' = 1, ..,N C = 3 



(28) 



The Heaviside smearing denoted by D Nv is particularly 
suitable since it cuts away the terms for k > N v , where 
the number of eigenvectors N v kept in the sum depends 
on the chosen width o" s =n,m,t<;- This choice of smearing 
reduces the number of needed inversions (per time slice, 
Dirac index and configuration) from the prohibitively 



1 A is a N X N matrix with eigenvalues A^ fc ' and eigenvectors , 
At;( fe ) = AWtjM (k = 1, ..,7V). 



FIG. 2. The spatial distribution *(r) $SU$ of the distillation 
operator □ * (|28|l constructed from the eigenvectors corre- 
sponding to the N v lowest eigenvalues of the Laplace opera- 
tor. The values are computed on each time slice of 49 config- 
urations at distances along the main axes and diagonals and 
plotted only until their respective symmetry points. Circles 
(black) , triangles (red) and stars (green) denote wide, middle 
and narrow sources (N v = 32,64,96), respectively. 



large number N^N C (needed for the conventional all-to- 
all approach) to a manageable number N v ~ O(100). 

Different truncations correspond to different effective 
smearing widths. We choose three smearing widths for 
quarks 

N v = 96 for s = n (narrow) , 
N v = 64 for s — m (middle) , 
N v — 32 for s — w (wide) , 



(29) 



which lead to the spatial distributions [20| of 



(30) 



shown in Fig. [2] 

We build each interpolator (|2"Tj) from quarks of the 
same width for all three widths. This enlarges the vari- 
ational basis and increases the possibility for optimal 
eigensets. 



F. Evaluation of the correlators 

The interpolators Ci_5 given in (|21[) arc linear combi- 
nations of quark-antiquark currents, which can be gener- 
ally written as 

gf c '(x',i) T a , a J*, c x (t,p) qf c (x,t), q,Q = u,d, 

(31) 

where the shape function Tit, p) incorporates the mo- 
mentum projection to p and the effect of covariant deriva- 
tives. Shape functions T for our interpolators (|21[) arc 
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given in (|A3[) of Appendix |XJ The pion-pion interpolator 
Og is a linear combination of products of two currents 

§U). 

After inserting the expression for smeared quarks q s 
of (|28p into interpolators (|2ip . all the contractions for 
C(tf,ti) can be expressed in terms of three quantities T, 
4> and r, analogous to the original proposal [20I ] which 
considered only one smearing width. Correlators are ex- 
pressed in terms of: 

• Dirac matrices T of size Nd x Nd- 

• The interpolator shape matrices (j)(t, !F) are square 
matrices of size N v x N v for an interpolator with a 
given smearing width N v 

cb k ' k (t,F) = J2 v&W^faPiv&Ht)' (32) 

x',X,c',C 

Our is related to $ in [2(| as $^ = fc ' fc r Q , Q . 

• The so-called perambulator matrices r k k (t' ,t) [20l | 
denote the propagators from source of shape v k (t) 
to the sink of shape v k (t') 

x' ,x,c ; ,c 

(33) 

Our correlators depend on the perambulators 
r(tf,ti), r(ti,tf), T(ti,ti), r(tf,t f ). These are 
in general rectangular matrices of sizes NdN? x 
N d Ni, N d Ni x N d Nj, N d N; x N d N* and N d Nj x 
NdNf respectively, where Nl^ = 32, 64, 96 denote 
the smearing widths of the source or sink. 

The analytic expressions for the needed contractions 
(Fig. [1]) in terms of T, r and <j> are given in Appendix El 

We precalculated and stored the perambulators 
T(tf,ti) from all source times slices U = 1,..,Nt = 32 
to all sink time slices tf = 1, ..,Nt = 32. This allows us 
to compute all needed contractions for C{tf,ti) straight- 
forwardly. We surr0 C(tf, U) over all initial time slices t, 
to decrease the relative errors on the resulting correlators 
C(t = t f -ti). 

We also sum over the results for the three p 
polarizations A = (0,0,1), (0,1,0), (1,0,0) for 
d = (0,0,0), or sum over the directions d = 
(0,0,1), (0,1,0), (1,0,0) for |d| = 1, and over the di- 
rections d = (1, 1, 0), (0,1,1), (1,0,1) for |d| = \/2- So, 
our final correlation matrices are 

C jk {t = t f -U)= E E C 3k{tf,U). (34) 

ti = l,..,JV T A or d 

These correlation functions finally enter the variational 
analysis ([T5)l to provide the energy levels. 



2 The sum plays the role of an average here. 



G. Finite N T effects and the "P+A" trick 

Our dynamical quarks have antipcriodic boundary 
conditions in time. Using the valence quarks with the 
same antiperiodic boundary condition in time, we find 
that the finite time extent N T = 32 (T = 3.96 fin) 
severely affects the eigenvalues \{t) near t ~ Nt/2 = 16. 
There are two major sources for this: 

• The 7r(pi)7r(p2) state receives contributions from 
both pions traveling forward or both traveling back- 
ward in time. But it also receives the contribution 
from 7r(pi) traveling forward and 7r(p2) traveling 
backward in time, and vice versa j35l |36|. As a re- 
sult, the cosh- type effective mass for some of the 
eigenvalues is not flat at t > 11. 

• In the pion channel, the ground state starts to dom- 
inate the second largest eigenvalue (and vice versa) 
at some moderate t [13, H3 • 

We use a previously applied trick, which effectively ex- 
tends the time direction to 2Nt = 64 by combining the 
periodic propagator Mp 1 and antiperiodic propagator 
M^ 1 (see for example J3(| UK]). All results in this paper 
have been obtained using the so-called "P+A" propaga- 
tors 

M -i ( , t s_[\[Mp\tf,U) + M2\t Jt U)] t f >U, 
M P+A t f ,t t) - | i [M -i {tfA) _ M -i {tfM tf < tt . 

(35) 

All our eigenvalues obtained from Mp^_ A agree with those 
obtained from A'/J 1 at t < 11. In the case of M^ 1 , the 
finite T effects seriously affect some of the eigenvalues for 
t > 11. In the case of Mp , A , the finite T effects never 
show up in any of the p eigenvalues for t < 16, which 
allows us stable fit ranges at least until t = 16. 

The "P+A trick" is not a valid field theoretic pre- 
scription, since the valence quarks do not have the same 
periodicity as the dynamical quarks (which remain an- 
tipcriodic in time). In practice, the pion correlators 
with zero momentum, for example, are perfectly consis- 
tent with periodicity 2Nt, i-e. they are proportional to 
e -m„t +e -mA2T-t) and keep f a u ing unti i t = 32 . We note 

that some of the nonzero momentum p correlators do not 
keep falling until i = 32, as would have been expected in 
the case of the proper field theoretic prescription. How- 
ever, none of the p correlators shows finite T effects for 
t < 16, which is the time window used for our analysis. 

III. COMPUTATIONS 

For the calculations presented here we use configura- 
tions generated for the study of reweighting techniques in 
the p-regime of chiral perturbation theory. A description 
of the normalized hypcrcubic smearing (nHYP smearing) 
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Nj x N T 

16 3 x 32 0.1283 7.1 0.1239(13) 1. 



k /3 a[fm] L[fm] ^configs m,r[MeV] 

280 266(3)(3) 



TABLE I. Configurations used for the current study. Nl and 
Nt denote the number of lattice points in spatial and time 
directions. For the determination of the lattice spacing a 
please refer to Sect. Mil The first error on m n is statistical 
while the second error is from the determination of the lattice 
scale. 



used in the dynamic fermion action can be found in [39| . 
Results from simulations with this action have previously 
been published in [29|, [3(| and the authors kindly pro- 
vided the gauge configurations used in this study. The 
action used to generate the gauge configurations contain- 
ing nf = 2 flavors of mass-degenerate light quarks is 
a tree-level improved Wilson-Clover action with gauge 
links smeared using one level of nHYP smearing. Ta- 
ble U lists the parameters used for the simulation along 
with the number of (approximately independent) gauge 
configurations used and the pion mass resulting from the 
determination of the lattice scale detailed in the next 
subsection. 

The gauge field obeys periodic boundary condition in 
time, while dynamical quarks are antiperiodic in time. As 
discussed in Sec. Ill Gl we compute and combine valence 
quark propagators with both antiperiodic and periodic 
boundary conditions. 

On each gauge configuration we calculate the lowest 
96 eigenvectors of the lattice Laplacian on every time 
slice using a standard 3-point stencil. Throughout, the 
gauge links are four dimensional nHYP smeared with the 
same parameters used for generating the gauge configura- 
tions: (01,02,03) = (0.75,0.6,0.3). For the calculation 
of the eigenmodes and the interpolating fields contain- 
ing covariant derivatives, we also experimented with ad- 
ditional three-dimensional link-smearing (using regular 
HYP smearing) and found only mild effects on the qual- 
ity of simple meson two-point correlators. We therefore 
opted to use no additional link-smearing. For the cal- 
culation of the eigenmodes we use the PRIMME pack- 
age In particular the routine JDQMR_ETol results 
in a fast determination for a small to moderate number 
(O(10) to 0(100)) of eigenmodes. For a larger numbcr 
of eigenmodes the Arnoldi/Lanczos method (and vari- 
ants) eventually outperform this method. For the meth- 
ods implemented in PRIMME we also tried a precon- 
ditioner using Chebychev polynomials, very similar to 
the method described in [18j . While this greatly im- 
proved the performance of some methods, our preferred 
method was largely unaffected and still outperformed 
all other PRIMME-methods for a moderate number of 
eigenmodes. 

For the determination of the quark propagators we 
use the df l_sap_gcr algorithm provided in Liischer's 
DDHMC package [4j], [42[ . Due to the large number 
of sources necessary for the distillation approach, an in- 



verter employing low-mode deflation techniques is espe- 
cially well suited. For the case presented here we ob- 
served a speedup factor of approximately five compared 
to a BiCGStab algorithm without low-mode deflation, 
while the computing time needed to generate the de- 
flation subspace was negligible compared to the actual 
calculation of quark propagators. Notice that this differ- 
ence gets more pronounced for the lighter quark masses 
needed for future studies at or close to the physical point. 

Statistical errors are determined with a single elimina- 
tion jackknife procedure throughout. When extracting 
energy levels we properly account for correlation in Eu- 
clidean time t by estimating the full covariance matrix in 
the given fit interval. For the covariance matrix we use 
a jackknife estimate which is calculated on the ensemble 
average onlj[f|. 

We determine the lattice spacing using the Sommer pa- 
rameter [44| ■ We extract the static potential from planar 
Wilson loops W(r, t) obtained on gauge configurations 
smeared with hypercubic blocking [45[ with standard pa- 
rameter values (0:1,0:2,03) = (0.75,0.6,0.3). The poten- 
tial is computed for each value of r from linear fits to 
log W(r, t) in the range t = 4 . . . 7 and then fitted to the 
lattice corrected form 



V(r) = A-\ \-ar + C 

r 



(36) 



in the range 1 < r < 7 or to the continuum form (i.e., 
C = 0) in the range 2 < r < 7. Both values agree within 
less than one standard deviation. The lattice corrections 
involves the lattice Coulomb pot ential [1/r] corrected for 
the hypercubic blocking fdll l47l| . To convert our num- 
bers to physical units (cf., Table HJ) we assume for the 
Sommer parameter the value ro = 0.48 fm and obtain 
a = 0.1239(13) fm. 



IV. RESULTS 
A. Pion results 

The pion energies are extracted from the variational 
analysis of the 6 interpolators given in (|2"4")) . The ex- 
tracted pion mass and pion energies for the two lowest 
nonzero momenta are given in Table [Hi along with the 
analytic predictions from the continuum and lattice dis- 
persion relations. 

B. Rho meson results 

1. Energy levels 

An example of the resulting correlators for interpola- 
tors 7T7r = Oq and qq = Oi, and their cross-correlators, 



This procedure has been referred to as jackknife reuse in [43 
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^1,2^1,2^1,2 


8-14 


1.57/5 


m n a =0.1673(16) 






(0,0,1) 


3 




12-17 


0.98/4 


0.4374(64) 


0.4268(65) 


0.4215(65) 


(1,1,0) 


4 




8-13 


1.31/4 


0.5823(46) 


0.5800(48) 


0.5690(47) 



TABLE II. The ground state pion energy extracted for three momenta: E is extracted from the variational analysis using the 
chosen interpolator sets, while E d,r ' are obtained using the ground state pion mass and the continuum and lattice dispersion 
relations (|10|l . 
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FIG. 3. An example of correlators for interpolators Ob = 7T7t 
and 0™ and their cross-correlators. 



are given in Fig. [31 

Given our 16 x 16 correlation matrices (j2"T]) . we ex- 
tracted the two lowest energy levels for a number of dif- 
ferent submatrices (i.e., interpolator sets) of dimension 
6 x 6 or less. Resulting levels for eight different choices 
of interpolator sets are shown in Fig. @] The extracted 
ground state energy is robust with regard to the choice 
of the interpolator set, while the first excited energy is 
robust only if the interpolator set includes the 7T7t inter- 
polator and if the correlation matrix is larger than 2x2. 
The first five choices include 7T7t in the interpolator ba- 
sis, while the last three do not. The first excited energy 
for d = (0,0,0) and d = (0,0, 1) has much larger errors 
and is often substantially higher if 7T7t is not in the set. 
On the other hand, it seems that the first excited energy 
in the case d = (1,1,0) can be extracted also without 
7T7r interpolator in the set. The choice set — 5 shows 
the result from the two-dimensional basis 7T7t = Oq and 
QQ = Ot, which was used by some previous simulations 
[lOL fl2l4l4l | . Figure 2] indicates that such a choice gives 
a reasonable estimate for the first excited energy in the 
cases d = (0, 0, 0) and d = (0, 0, 1), while it gives a much 
higher energy for the first excited state with d = (1, 1, 0). 
Our study shows that a basis larger than 2 x 2 is needed 
to extract the first excited level in this case. 

Given that our lowest two energy levels are robust with 
respect to the choice of interpolator set provided the set is 
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FIG. 4. The lowest two energy levels (circles denoting the 
ground state, squares the 1st excited state) extracted using 
different submatrices (interpolators sets) of the full 16 x 16 
correlation matrix (|21[) . all for to = 4. The horizontal dashed 
lines indicate the energy values for two noninteracting pi- 
ons. The interpolators Ois have qq valence structure, while 
Ofj = mr. All interpolators in this plot are composed of nar- 
row quarks q s = q n , with the exception of interpolator set 3 
which is O^O^O^Ol. In order to make different interpola- 
tor choices comparable, we use the same fit range t = 7 — 10 
in the one-exponential correlated fit for the purpose of this 
figure (with exception of E2(d = 0) obtained for t = 5 — 7). 



large enough and contains the 7T7t interpolator, we present 
the final interpolator set choices in Table Hill The corre- 
sponding effective masses for our preferred interpolator 
choices are shown in Fig. [SJ The final values for the six 
energy levels in Table Mil are extracted using correlated 
two-exponential fits with to as indicated in the table and 
starting at a rather small time separation t. We veri- 
fied that the extracted levels agree with results obtained 
from one-exponential fits starting at larger t and using 
to = [2,5]. 



2. Phase-shifts and resonance parameters 

Each of the energy levels of Table IIIII gives the value 
of the scattering phase shift S(s) at a different pion CMF 
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level n 


to 


Interpol. 


fit range 


E n a 


X 2 /d.o.f. 


ap* 


s a 2 


8 


(0,0,0) 


1 


2 


01,2,3,4,6 


3-18 


0.5107(40) 


6.10/12 


0.1940(29) 


0.2608(41) 


130.56(1.37) 


(0,0,0) 


2 


2 


^1,2,3,4,6 


3-12 


0.9002(101) 


0.85/6 


0.4251(58) 


0.8103(182) 


146.03 (6.58) [*] 


(0,0,1) 


1 


2 


Cl,2,3,4,6 


4-16 


0.5517(26) 


4.06/9 


0.1076(36) 


0.1579(29) 


3.06 (0.06) 


(0,0,1) 


2 


2 


Cl,2, 3, 4, 6 


4-15 


0.6845(49) 


3.10/8 


0.2329 (40) 


0.3260(69) 


156.41(1.56) 


(1,1,0) 


1 


3 


C™2,3,6 


4-12 


0.6933(33) 


4.33/5 


0.1426(42) 


0.1926(49) 


6.87(0.38) 


(1,1,0) 


2 


3 


Cl,2,3,6 


4-12 


0.7868(116) 


2.34/5 


0.2392(101) 


0.3375(191) 


164.25(3.53) 



TABLE III. Final results for the lowest two p energy levels, all obtained using 2-exp correlated fits with given % 2 /d.o.f.. The 
choice of interpolator basis (|2ip is indicated. The pion momenta ap* in the CMF and scattering phases 5 are obtained using 
the lattice dispersion relation (|10[) and a in Table HT1 The state Ez(P = 0) is above the 4n threshold and is denoted by a 
star. 



momentum p* . We employed the lattice dispersion rela- 
tion (fl"U| to get p* = ^-q and used the phase-shift for- 
mulas in Sect. Ill Al to get S(q 2 ). Our results, including 
jack-knife error estimates, are also given in Table [TTTT 

The resulting phase-shift is related to the relativistic 
Breit-Wigner form for the elastic p-wave amplitude in 
the resonance region [l| 



ai = 



m 2 + iy/s r(s 



(37) 



where s = Eq M is the Mandclstam variable and m 2 is 
the resonance position. Relation (|3"T|) can be conveniently 
written for later use as 



f~s r(s) cot 8{s) = m p — s 



(38) 



and the decay width T(s) is expressed in terms of the 
coupling constant g p7r7r , taking into account the tttt phase 
space [S Hi] 



r(*) 



P* 3 9%* 
s 6ir 



(39) 



where the p width T p = T(m 2 p ) is evaluated at the reso- 
nance position. 

The final relation, the so-called effective range formula, 
combines (|38I39[) and is valid in the elastic region s < 
(4m,,) 2 , 



P 



cot S(s) = -^-(mj - s) 

9 P7T7T 



(40) 



It allows a linear fit for the two unknown parameters 
Gir/g 2 ^ and 6tt m 2 / 'g 2 w ^. Values of s, p* and 5 for the 
energy levels E n are given in Table IIIII and appropriate 
combinations are plotted in Fig. [5] In the fit and in 
the figures we do not include the first excited state with 
P = 0, since this lies above the Air inelastic threshold. 

Figure [6] shows the result of the linear fit to the data, 
giving our final result for gp^^ and the mass of the p 
resonance (at our = 266(3)(3) MeV), 



g p ^ = 5.13(20) , 
m p a = 0.4972(42) 



(41) 



Figure [7] exhibits the corresponding phase-shift in the 
resonance region. The values (|4ip are obtained using the 
lattice dispersions relation (fTUj) . Given the systematic 
uncertainty with simulations on a single ensemble, they 
agree reasonably well with the results <7 p7r7r = 5.60(18) 
and m p a = 0.4833(41) obtained using the naive disper- 
sion relation. 

The value of the coupling (|4"Tj) is near the experimental 
value gp^v « 5.97. Our coupling is also compatible with 
the results in [l(J [l2[ within the errors given there. Note 
that fToL fl^] computed the coupling at a larger pion mass. 
In [13j a larger value (7 p7r7r = 6.77(67) and a substantially 
larger m p = 980 MeV are observed at similar pion mass 
m-x = 290 MeV. Our m p is close to the prediction of the 
unitarized one- loop@ Chiral Perturbation theory (ChPT), 
which leads to about m p ~ 800MeV at m w ~ 266MeV 
[5(l[5l|. We also compared our S(s) with the prediction 
of unitarized ChPT, recalculated for our = 266McV 
by the authors of [52]: we find good agreement for y/s < 
Trip and reasonable agreement with one-loop results for 
yfs > nip. 

Since the width is crucially influenced by the 7T7r-phase 
space, this number derived for our pion mass comes out 
significantly smaller than the experimental value, so we 
present only gp-x-n- This dimensionless coupling is ex- 
pected to be almost independent of m w [5l| , which was 
also explicitly verified in a study for several pion masses 

M. 



V. CONCLUSIONS AND OUTLOOK 

Extracting scattering phase-shifts and resonance prop- 
erties is one of the most challenging problems in hadron 
spectroscopy based on lattice QCD. We combine several 
sophisticated tools to approach this problem: Liischer's 



m p = 792(7) (8) MeV 



4 The two loop result strongly depends on a number of poorly 
known Low Energy Constants, which are fixed in l5ll also by 
using the lattice data on m p , so the comparison to the two-loop 
result is not appropriate. 
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FIG. 6. Our data for ((ap*) 3 /^, 



cot 5(s) as a function of 



* , fitted to straight line behavior according to (|40|1 . The fit 
has x 2 /d.o.f.=7.42/3 and gives g p7r n = 5.13(20) and m p a — 
0.4972(42). The states E n (d) corresponding to various points 
can be deduced by the value of s in Table [TTT1 The plot data 
is shown in units of the lattice spacing. 
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FIG. 5. The effective energies observed in the three momen- 
tum frames (0,0,0), (0,0,1) and (1,1,0), based on diagonaliza- 
tion of a correlation matrix with 4 or 5 operators, listed in 
Table HTT1 The horizontal bands indicate the resulting energy 
levels derived from two-exponential fits to \i(t) as discussed 
in the text. The dashed lines indicate the noninteracting two- 
pion levels as determined from the energies Ef a l'a in Table fTTl 



phase-shift relations for finite-volume lattices, moving 
frames and variational analysis of correlation matrices, 
where a number of quark-antiquark and 7T7t interpolators 
with quantum numbers I(J PC ) = 1(1") are used. All 
needed contractions are evaluated using the distillation 
method with the Laplacian Heaviside smearing of quarks. 
We find that these tools lead to precise values of the p- 
wave phase-shift for 7T7t scattering at five values of pion 



FIG. 7. The p-wave phase-shift values compared with the 
result from the fit to ()40[) in Fig. [5] for gp-mr = 5.13 and 
m p a = 0.4972. The states E n (d) corresponding to various 
points can be deduced by the value of s in Table IHII 



relative momenta in the vicinity of the resonance. This 
allows a determination of the p resonance parameters m p 
and T p at our value of m n . 

The simulation is performed on an ensemble [29|, [3(| of 
280 gauge configurations with two mass-degenerate dy- 
namical clover-improved Wilson fcrmions. The pion mass 
TOtt is roug hly 266 McV, the lattice volume V is 16 3 x 32 
and the spatial extent of the lattice is L ~ 1.98 fm. The 
exponentially suppressed finite-volume corrections may 
not be completely negligible at our m^L ~ 2.68 and fu- 
ture simulations will have to improve on this. Larger lat- 
tices will necessitate stochastic estimation techniques to 
avoid the unsatisfactory scaling of full distillation with 
the lattice volume. Such a method has recently been 
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provided in (l8|. In the present study we calculated the 
quark propagation by calculating the distillation peram- 
bulators on all time slices, which is not very economical 
and only feasible in small volumes. 

Along the way, we explore how well the lowest two en- 
ergy levels can be obtained without the tttt interpolators 
in the variational basis. We also propose how to treat 
interpolators of different smearing widths in the same 
variational basis within the distillation method. 

We demonstrate that a relatively accurate determi- 
nation of the resonance parameters is possible with 
present day techniques, within the limitation of small 
m^L. For our pion mass we obtain the resonance 
mass m p = 792(7) (8) McV and the p — >• tttt coupling 
9pTTK = 5.13(20), which is close to the experimental value 
g^jjP w 5.97. We prefer to give the coupling, since the 
actual width T p is strongly affected by the phase space, 
which is small due to the large value of our pion mass. 

Following the pion, the rho is the most prominent me- 
son. With sharpened tools it is now becoming possible to 
analyze its decay properties. The present study of the p 
resonance gives us confidence that similar techniques can 
be applied to also extract the resonance parameters of 
some other hadronic resonances and we intend to pursue 
research along these lines in the near future. 
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Appendix A: Contractions in the distillation method 

Here we provide the analytic expressions for correla- 
tors C(tf,ti) = (Of(tf) Oi(U)) that follow from general 
quark antiquark interpolators with \I, 1 3} = |1,0) (with 
examples given by O1-5 in (|2ip) 



Of (if) = 77|[M*/) r ° P) «-,(*/) - d }} > 

of{u) = -L [u Sz (t t ) r° r°(u, -p) u St (u) - {u d}] , 

(Al) 

and general meson-meson (MM) interpolators with 
\I, I 3 ) = |1, 0) (with example given by O e in {2TJ) 



O 



1 



MM 

u Sf (t f ) T 2 f J"/(t/,p 2 f) d Sf (t) - {u^ d}] , 

x 



1 



of M (u) = ^=[u Si {ti) r M ^(ti.-pu) d Si (t, 

d Si (ti) r 2l J r t (< l ,-P2i) u Si (t) - {uo d}] . 

(A2) 

The subscripts Sj, Sf = n,m,w denote the smearing 
width of the sink and source. The superscript "0" de- 
notes that r° and T° apply to qq interpolators, while T 
and T without superscript apply to meson-meson inter- 
polators. 

The shape functions J- (|31l) for our interpolators (f2"Tj) 
are of three types 

no V: J^(*,p) =<WW ipx , 

for V : 7&(t,p) = \ [e^^Ox'xW - ^ it x(i)c ipx ] , 
forVV:^(i,p)= Y. (^)x'xo(^ ipXo (^)x°oxW, 



i=l,2,3 



(A3) 



and we use the first choice (without V) within our meson- 
meson interpolators. 

The contractions in Fig. [1] are expressed in terms of 
the perambulators t (|33[) . interpolator shape matrices 4> 
((32)) and Dirac matrices T, which are presented in Section 
IIIFI of the main text. The analytic expressions for the 
contractions are 



(Of(t f )Or(U)) = C Fi °- lb (t f ,U) = -Tr[ r(U,t f ) T° f 0(t/,^(P)) r(t f ,U) T» <f>(U,F?(-P)) 



(A4) 



Tr[ r(U,t f ) T lf <t>(t f ,F(pit)) r(t ft t f ) T 2f J>(t f ,F(p 2f )) r{t f ,U) T° <j>(t u 7?(-P)) ] 

- {pif o P2f , Ti/ o r 2/ } . 



(A5) 



- Tr[ T(t f ,ti) r u 4>(ti,F(-Pii)) r(U,U) T 2l <f>(U,F(-p 2i )) r(fc,f/) T° f <j>(t f ,F° f (P)) ] 

+ {pii <-» P2i , Tu <-> r 2 j . 

{0 MM {tf)0 MM {ti)) = G^itf^ + C^it^U) . 

C^ Aa (tf,U) = Tr[ T(ti,t f ) T lf fttf^fipu)) r(t f ,U) T u ^(-pii)) ] 
x Tr[ T(U,t f ) T 2/ <t>{t f , J"/(p 2 f)) T(t/,t<) T 2l Ji(-pai)) ] 

- {pii <-» P2i , r H o r 2i } , 

C2 9 - lc (t f ,U) = Tr[ T (U,t f ) T 2/ 0(t/,.F/(p a ,)) T(t f ,t f ) T lf c/>(tf, J-/( P if)) 
• T(t/,t<) r 2l 0(«j,7i(-pai)) rfotj) r u ^(-pn)) ] 
+ Tr[ T(U,t f ) T lf <j>(t f ,F f (p 1{ )) T(t f ,t f ) T 2/ <t>{t } , Fffau)) 

■ r(t f ,ti) T u 0(ti,J"i(-pii)) rfc.ii) T 2l <j>(U, Ji(-pai)) ] 

- {ph <->■ P2i , r H «-» r 2 j . 



and can be generally used for the interpolators of the 
form (fXI) or (|A2[) , or their cross-correlators. 
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